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ABSTRACT 

The velocity, position, and action variable evolution of a tidal stream drawn out of a star cluster in a 
triaxial isochrone potential containing a sub-halo population reproduces many of the orbital effects of 
more general cosmological halos but allows easy calculation of orbital actions. We employ a spherical 
shell code which we show accurately reproduces the results of a tree gravity code for a collisionless star 
cluster. Streams from clusters on high eccentricity orbits, e > 0.6, can spread out so much that the 
amount of material at high enough surface density to stand out on the sky may be only a few percent 
of the stream’s total mass. Low eccentricity streams remain more spatially coherent, but sub-halos 
both broaden the stream and displace the centerline with details depending on the orbits allowed 
within the potential. Overall, the majority of stream particles have changes in their total actions of 
only 1-2%, leaving the mean stream relatively undisturbed. A halo with 1% of the mass in sub-halos 
typically spreads the velocity distribution about a factor of two wider than would be expected for a 
smooth halo. Strong density variations, “gaps”, along with mean velocity offsets, are clearly detected 
in low eccentricity streams for even a 0.2% sub-halo mass fraction. Around one hundred velocity 
measurements per kiloparsec of stream will enable tests for the presence of a local sub-halo density 
as small as 0.2-0.5% of the local mass density, with about 1% predicted for 30 kiloparsec orbital radii 
streams. 

Subject headings: dark matter; Galaxy: structure; Galaxy: kinematics and dynamics 


1. INTRODUCTION 

Tidal star streams pulled out from star clusters have 
been extensively studied in simulations where the galaxy 
is modeled as an analytic potential, often taken to be 
spherically or axially symmetric for simplicity (|McGlynn 


1990; Joh nston et al.|1 996; Dehnen et al.|200 4; Kiipperet 
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& Burkert| |2015|) . These simnla- 


al. 120081 |Varghese et'al.|2011||Kupper et aL|2012||Amor- 


2014; Fardal et ah| 20141 |Ngan &; Carlberg 2014 


tions have described how tidal heating ejects stars into a 
stream. There is a good analytic understanding of how 
sub-halos interact with thin streams in angular momen¬ 
tum conserving potentials to create ve locity pertu rba 
tions which open up to gaps in a stream (Yoon, Johnston 
& Hogg||201l| |Carlberg| 2012| |Erkal & Belokurov||2015p. 
Simulations that begin with a cosmological halo show 
considerab ly more orbit diffusion than smooth analytic 
potentials flBo naca et al.|2014 Pearson et al.|2015 Ngan 
et al. 2015a bj which is not yet well understood. How- 
ever, the increased orbital freedom is certainly related to 
the rich orbit str uctu re in static triaxi al potentials (de 


Zeeuw 1985 Binney & Tremaine 2008) which is bound 
to"play a role in the aspherical and non-stationary cos¬ 
mological galactic halos in which tidal streams develop. 

In principle, kinematic measurements along a stream 
can tightly constrain the full phase space parameters 
of the highly redundant orbits of the stream particles 
and ca n serve as a precisio n probe of a galactic po - 


tential (Johnston et al. 1999 
The initial velocities of streams 


Binney||2008| |Bovy||2014|). 
is of stars that are tidally 
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stripped from star clusters are calculated for a given orbit 
and potential with collisionless n-body codes. The stars 
emerge with highly correlated orbital velocity substruc¬ 
ture (Carlberg 2015). The stream width scales with the 
tidal radius and the length of scales with th e tidal radius 
and the age of the stream (Johnstonll 19981). In a spher¬ 


ical potential angular momentum and radial action are 
always conserved, but tidal mass loss features blur out 
with orbital overlap down the stream. However as stars 
move away from the progenitor on slightly different or¬ 
bits in which the gravitational field has variations due to 
aspherical structure in the halo the streams can disperse 



et al. 

2014) 

Ngan et al. 2015b 

). Relative to a spherical 

poten 

tial the streams are almost everywhere increased in 


width, with some portions greatly blurred. However, a 
few se gments of high or e ven enhanced density can re¬ 
main (Ngan et al. 2015b). The view of the stream in 
the dynamical action variables, which are adiabatically 
conserved quantities of a potential, can provide insight 
i nto what is causing the stream to change its appearance 
(jlbata et al.||2002|. 

The purpose of this paper is to quantitatively analyze 
the dynamical behavior of tidal streams in moderately 
aspherical potentials that also contain a sub-halo pop¬ 
ulation. The analysis will be done within a simplified 
potential which is able to reproduce many of the fea¬ 
tures of a cosmological halo but allows the calculation 
of action-angle variables at the beginning and end of the 
simulation. A simplified n-body code that reproduces 
the tidal stream from a standard n-body code is used 
to speed up the simulations. The changes in angular 
momentum and the radial action along a stream as a 
function of the halo shape and sub-halo content appear 
as changes to the stream in position and velocity space. 
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These results will be useful to advance the understand¬ 
ing of the mass structure and contents of galactic halos 
as more and better kinematic data become available and 
serve as a guide to more realistic simulations. 


2. ORBITS IN TRIAXIAL ISOCHRONE 
POTENTIALS 

We want to analyze the orbit changes of a stellar 
stream in a galactic halo with sub-halos. Action vari¬ 
ables, as conserved quantities, provide direct insight into 
what has changed. Henon’s isochrone p otential has sim¬ 
ple algebraic expression for the actions (Henon 1959ajb ) 
and provides a reasonable match to a galactic rotation 
curve over the range of radii where we want to study 
stream dynamics. One fairly standard way to create a 
triaxial potential is to apply a linear scaling to the Carte¬ 
sian coordinates of a spherical potential to create a tri¬ 
axial potential. The simplicity of the axis scaling comes 
at the cost of losing the analytic properties of the orbits 
in the spherical potential. A number of axisymmetric 
potential s with good analytic properties hav e been de¬ 
veloped (|Eyre & Binney||2011| |Binney||2014|) but these 
still do not have the triaxiality, or more generally, the 
asphericity, of cosmological halos. 

With the goal of allowing the rich orbital structure of a 
triaxial potential yet retaining some of the analytic prop¬ 
erties of the isochrone, we introduce a time dependent po¬ 
tential, which for cosmological halos should be a reason¬ 
able approximation. The potential starts as a spherical 
isochrone, makes a smooth transition to a triaxial form, 
then remains as a triaxial for a relatively long interval, 
then reverses to become spherical again. The scale ra¬ 
dius and the mass remain constant, the shape being only 
thing that changes. We will leave the x axis unaltered, 
but compress the y — z axes, which alters the density 
profile, increasing the density in the central region and 
generally leading to negative densities at large radius on 
the most compressed axis. We will also smoothly turn 
on and off a small fraction of the mass in an identically 
distributed sub-halo population. 

The potential’s spherical initial and final form is a de¬ 
vice to allow the analytic actions to be calculated and is 
not intended to be realistic. The action variables of the 
spherical isochrone will be valid quantities to calculate 
at the beginning and end of the simulation, although 
they are not conserved quantities over the simulation. 
We expect that there will be smooth changes in the ac¬ 
tion variables associated with the turn-on and turn-off 
of the triaxiality, and local changes as a result of sub¬ 
halo encounters, the two being sufficiently different that 
it might be possible to untangle the effects. Based on 
these considerations the triaxial isochrone potential used 
here is, 


$(x,y,z,t) = 


-GM{t) 


b+ fb 2 ~+~x 2 ~+y\y 2 (t) + Tq z 2 (i) 


• ( 1 ) 


The q y (t) and q z (t ) are dimensionless functions which 
begin and end at unity, at which time the potential is 
spherically symmetric. We largely focus on compres¬ 
sions where the q values are less than one, but cosmolog¬ 
ical halos also undergo small expansions during mergers 
with other halos as they. Our standard compression has 


q y = 0.9 and q z = 0.95 which leads to small negative 
densities on the y axis beyond y = 6.51, with a negative 
peak at y = 8.83 with about 10 -4 of the central density. 
The negative densities occur outside of the orbits that 
will be studied. The gravitational force is always inward 
directed. 


2.1. Time Variation of Potential Shape 

The potential asphericity turn-on takes place over a 
time £ r , after which it remains on for a time £/, using 
the arbritrary function, 

A(t) =2(t/t r ) 2 0 < t < t r / 2, 

=1 - 2(1 - t/t r ) 2 t r /2 <t<t r , 

=1 t r < t < tf + t r . (2) 

For t > tf + t r the rising function is reversed and A(t) 
returns smoothly back to zero. The rise time, £ r , is cho¬ 
sen to reasonably slow, typically about 3 radial orbits, 
which is neither in the impulsive regime nor truly in the 
adiabatic regime. The time spent as a constant triaxial 
potential, £/, is typically about 15 radial orbits, com¬ 
parable to the time over which thin streams have been 
developing in the Milky Way halo. 

All calculations here use units with GM = 1 and 6=1. 
We usually start the calculation with a cluster placed at 
an x position of 4 units. We characterize the starting con¬ 
ditions with the angular momentum, L, relative to the 
angular momentum of a circular orbit in the spherical 
potential, L c . Orbits with L/L c of 0.4 and 0.7 have ra¬ 
dial orbital periods are 29.2 and 35.3, respectively, in the 
spherical potential when started at r = 4. For the time 
variation function A(t) we select t r = 80 and tf = 440, 
so that the potential becomes spherical again at time 600 
and the simulation terminates. 

The co-ordinate distortions are evaluated as q y (t ) = 
1 + a y A(t ) and q z (t ) = 1 + a z A(t). Positive a y ,a z cor¬ 
respond to an expansion of the potential, that is, re¬ 
duced gravitational attraction at locations off the x axis 
and negative values correspond to a compression. Both 
expansion and compression are possible in cosmological 
situations, although compression likely dominates and is 
relatively small, since accreted mass spends little time 
in the central region and largely deposits mass near the 
current virial radius of the halo. Cosmological halos also 
can have a rotation of the potential shape, whic h we will 
ignore at this time. Louis & Gerhard (1988) studied 
radial perturbations to the isochrone potential finding 
that higher order resonances quickly appeared along with 
complex orbits, leading us to expect behavior of that sort. 


2.2. Orbiting Sub-Halos 
Sub-halos are present in large numbers in the dark 


matter halos created in cosmological simulations ( 

Klypin 

et al. 19991 Moore et al. 

1999; piemand, Ku 

men & 

Madau 2007 Springel et al 

L 2008j Stadel et al. 

2009). 


hence has nearly a constant mass per decade of mass, 
m 2 dN/dm oc m 01 . The mass fraction of sub-halos varies 
steeply with radius relativ e to t h e local d ark matter den¬ 
sity in cosmological halos (Spr ingel et al.| |2008). For the 


radial range around 30 kpc where the few known thin 
streams have been found the models predict about 1% of 
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the mass is in sub-halos, with the sub-halo mass fraction 
dropping to about 0.1% at a nominal solar radius of 8 
kpc. 

For the interests of this paper the practical issue is to 
have a well defined sub-halo population over the radial 
range of the star cluster orbit. Conceptually our ap¬ 
proach is to convert a small fraction, f s h, of the mass of 
the fixed isochrone potential into a set of sub-halos with 
the same collective massive distribution as the isochrone. 
The sub-halos are free to orbit within the isochrone po¬ 
tential as test particles, which requires that they be ran¬ 
domly drawn from a distribution function which gener¬ 
ates the isochrone potential. Conveniently analytic forms 
are available for t he isochrone ’ s mass and en ergy distri¬ 
bution functions (|Henon||1960| |Binney||2014|) . Therefore 
on the average there is no change in the summed sub-halo 
plus diminished external potential, and, on the average, 
the global and local sub-halo fractions are equal every¬ 
where for the initial spherical potential. Although the 
density distribution of the isochrone extends to infinity, 
59% of the mass is within the r = 4 orbit where we start 
the star clusters. Hence, little computational effort is 
wasted on the relatively few sub-clusters at larger radii 
which do not have orbital pericenters within the range 
of radii that the star cluster explores. The sub-halos re¬ 
spond to the undiminished triaxial isochrone potential 
but not each other. The stars of the cluster and the tidal 
streams feel gravitational forces from all other stars, plus 
the isochrone reduced by (1 — f s hA(t )) plus the sum of 
the forces from all the sub-halos which vary as f s hA(t). 
where the A(t) factor, Equation [2] which rises from zero 
and returns to zero at the end ofthe simulation. 

All the sub-halos are given identical masses, thereby 
allowing us to examine mass dependent effects. That is, 
once a number, N s h , of sub-halos is chosen, the sub-halo 
all have mass, m s h = f s h/N s h • Each sub-halo is_mod- 
elecl as a Herquist potential,— Gm/(r ± a), (Hernquist 


1990), with a scale radius a s h = tym s h, as expected for 


cosmological sub-halos of constant scale density. The 
dependence of the results on the scale radius of thejsub- 


h alos are visible in the simp lified analytic analysis (Carl- 
berg| 2015 Erkal & Belokurov] 2015), where there is little 


change beyond the scale radius, but the induced veloc¬ 
ities increase for stars that have encounters that take 
them within the scale radius of a sub-halo. 


2.3. Test Particle Orbits Without Sub-halos 

Simple orbit integrations in the model potentials illus¬ 
trate some of the potential’s basic properties. We start 
a set of 10,000 particles with no self-gravity at x = 4, 
y = z = 0 with a uniform distribution in v y between zero 
and the circular value at the chosen starting radius. The 
value of x = 4 is chosen to help avoid the core region 
where there is little effective tidal force. The particle 
v x and v z velocities are given a small Gaussian random 
value, 0.002 units, about 0.5% of the circular velocity, to 
allow the particles to explore mildly off plane orbits. The 
orbits are integrated with a leapfrog integrator with step 
sizes of 0.004 units, which is nearly 10 4 steps per orbit, 
conserving the actions to better than one part in 10 6 . 

Orbits are evolved in four different versions of the triax¬ 
ial isochrone, with a y ,a z values of ±0.05 and ±0.1. The 
differences of the actions of the particles calculated at the 
start and end of the orbit integrations are plotted as a 



L 


Fig. 1.— The changes between beginning and end of the simu¬ 
lation in the particle’s total angular momentum, L (blue), radial 
action J r (red) and angular momentum about the axis perpendic¬ 
ular to the initial orbital plane, J z (green), as a function of the 
initial L in a potential with a y = —0.05 and a z = —0.10. The 
particles are started at r = 4 with a uniform angular momentum 
distribution between 0 and the circular value. All orbits are started 
in the plane with an initial vertical and radial velocity dispersion 
of 0.002 units. A point below the downward sloping line drawn 
from the origin would have had its orbital direction reversed. 
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Fig. 2.— Same as Fig. [I] but for a y = —0.10 and a z = —0.05. 


function of the initial orbital angular momentum in Fig¬ 
ures [TJ [2j [3] and [4] Figure [l] shows changes roughly as one 
might expect as a function of the orbital angular momen¬ 
tum (green points). That is, the range of initial orbital 
angular momentum means that particles have different 
orbital phasing relative to the changing axes. Some par¬ 
ticles experience a net positive torque and others where 
the torques lead to a net loss. At high L, hence low ec¬ 
centricity, the torques largely cancel out, leaving little 
net change. 

In contrast to the simple behavior of Figure 0 Fig- 
ures[2][3]and [4]show nearly chaotic b ehav ior over a range 
of orbits (Siegal-Gaskins & Valluri|2008|). The change in 
potential axis length remains the same as in Figure [l] 
but with the yz axes switched, Fig. [2] or, made pos- 
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Fig. 3.— Same as Fig. [T]but for a y = 0.05 and a z = 0.10. 



Fig. 4.— Same as Fig. [T]but for a y = 0.10 and a z = 0.05. 



L 


Fig. 6.— Orbit integrations for the potential of Fig. [2]with 1000 
sub-halos that total to 1% of the halo mass. The grey lines show the 
same simulation with the sub-halos reduced in mass by two. Here 
the particles are initiated with no additional random velocities. 
The inset shows the same simulation with no sub-halos. 


itive, hence allowed to expand, Figs. [ 3 ] and [4] The fig¬ 
ures show that nearby orbits have dramatic differences in 
the changes to their action variables for higher eccentric¬ 
ity orbits. Low eccentricity orbits at high angular mo¬ 
mentum are generally much less affected, although the 
vertical action can still be dramatically affected. This 
behavior is often seen for orbits that have higher order 
resonances between orbits and the rate of potential shape 
change. 

If the potential changes are made very slowly we ex¬ 
pect that the action changes for regular orbits will be- 


come adiabatic an d become exponentially small ( |Binney 
fe Tremaine|2008 ). In Figure [5] the A L and A J 7 


actions 

diminish as expected but the changes in the A J z values 
remain large. Overall, the rich range of orbital behavior 
of the triaxial isochrone makes it a simplified, but inter¬ 
esting, stand-in for a cosmological galactic halo potential 
that can capture realistic star stream behavior beyond a 
spherical halo. 



Fig. 5.— Same as Fig. [2j with an increased vertical axis, but 
for t r = 2400, which is 30 times slower than the standard rate of 
potential change. The green dots are the changes in J z . 


2.4. Test Particle Orbits With Sub-halos 

The effects of adding sub-halos to the simple orbit inte¬ 
grations of Figure |2| are shown in Figure [6j where 1% of 
the halo mass transfers to 1000 identical mass sub-halos. 
That is, each sub-halo has a mass 10 -5 of the host halo. 
The test particles are initiated with no additional ran¬ 
dom velocities to avoid the confusion of orbit diffusion. 
The outcome with no sub-halos is shown as the inset in 
Figure [6] 

The sub-halo interactions occur on the scale of the ra¬ 
dius of the sub-halos, which is approximately 0.02 units 
for the case studied. The sub-halos act to create fairly 
abrupt changes in the action variables, spread out over 
the orbits in accord with the density profile of the halo. 
Reducing the mass of each sub-halo by a factor of two, 
which reduces their scale radius by v^2, we find that the 
sharp changes in J z for L > 0.8 are one half the height 
of the higher mass sub-halo results, indicating that sin¬ 
gle sub-halo encounters dominate for test particle orbits 
in this range over the duration of the simulation. For 
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Fig. 7.— The mass in the cluster as a function of time on an 
L/L c = 0.4 orbit evolved with an m = 10 -7 with a full n-body 
code (solid line) and with a shell code (dashed line). 
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Fig. 8. — The angular momentum (blue) and radial action (red) 
relative to the cluster center, AL,AJ r , respectively of particles 
as a function of their ejection time for the cluster evolved in the 
Gadget2 n-body code. 


L <0.8 the particles encounter sub-halos multiple times 
leading to much larger J z changes. The highest eccentric¬ 
ity orbits, those with L < 0.3 are largely unperturbed by 
sub-halos, likely because of their high velocities through 
the central regions. Overall the amplitude of the changes 
in the actions due to sub-halos are comparable to the 
changes due to the triaxial potential. We can foresee 
that disentangling sub-halo effects from potential effects 
will not be straightforward as might be hoped. 


3. SIMULATION OF A SMALL STAR CLUSTER 

We create a small, self-gravitating satellite of mass of 
m < 1 (m ~ 1CV 8 , where the host halo h as Mj = 1) 
from a King model distribution function (King 1966). 
The cluster is created with its own dimensionless units 
of velocity and position, then scaled down in size to fit 
its zero density radius within the estimated tidal radius. 
The velocities are scaled to the desired mass and new size. 
Our standard star cluster model has a central potential, 
V’(O), relative to the central velocity dispersion, cr, of 
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Fig. 9.— The same model as Fig. [8] but evolved with the shell 
code. 

—^(0)/cr 2 = 4. Increasing —^(O)/cr 2 decreases the core 
radius relative to the radius where the model density 
goes to zero. The outer density profile has a slope nearly 
independent of the core radius, so the mass loss into the 
tidal stream has a weak dependence on the concentration 
of the model star cluster. Consequently the star cluster’s 
core radius only begins to have a significant effect on the 
tidal stream when mass loss is so large that the core is 
eroded. 

The outer radius of the star cluster model is scaled to 
be equal to the tidal radius, r[m/ (3Mj(r))] 1 / 3 , where r is 
the starting radius of the orbit. The small adjustment to 
the tidal radius fo r a distributed mass potential (Binney 


& Tremaine 2008) is not included, nor is any allowance 
for a non-circular orbit, which leads to some modest in¬ 
ternal changes of the star cluster model in response to 
the tidal field in its first orbit. 

The standard model uses 10 5 particles in the initial star 
cluster. The star particles in the model star clusters are 
subject to gravitational forces from other star particles, 
the host galaxy gravitational potential, and any orbiting 
sub-halos that may be present. The rate of mass ejection 
varies around the orbit and leads to complex, strongly 
correlated, velocity differences (or action variable o ffsets) 
between the stream and the star cluster (Carlberg|2015). 
The mass loss in the systems of interest is driven by tidal 
heating and is well modeled with a collisionless n-body 
code. As a reference, we use the gravitational and orbital 
evolution componen ts of the well known and character¬ 
ized Gadget2 code (Springel 2005). 


3.1. An n-body shell code 

The precise description of the evolution of a star cluster 
in the tidal field of a host galaxy remains an important 
research problem as a generalization of the dynamics of 
an isolated star cluster. Current work allows for a real¬ 
istic range of star masses, their mass evolution and their 
remnants, the presence of binaries and higher multiplic¬ 
ity systems. The evolution of collisional star clusters in 
tidal fields generally either use some form of gravitational 
calculation sp litting to handle the dynamics on different 
time scales (Fujii et al. 2007)_ or i nteractio ns via statis¬ 
tical approximations (Sollima & Mas trobuono Battisti 
2014 Brockamp et al. 2014]r Fukushige & Heggie (2000) 
demonstrate the complexity of star cluster orbits near 
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the tidal radius that exists for cluster in a circular orbit. 

Two well known thin star streams are Pal 5, which 
ema nates from an exceptiona lly low concentration clus¬ 
ter ( Odenkirchen et al. |2Q0l [) and GD-1, which has no 


known progenitor (Grillmair & Dionatos 2006). Given 
the complexity and high computationaFcosts of study¬ 
ing a collisional star cluster, and, the observational fact 
that the importance of star-star encounters in the pro¬ 
genitor clusters is either relatively low, or, unknown, it is 
currently conventional to study tidal streams as emanat- 
ing from tidally heated collisionless star clusters (Dehnen 
et al. 2004). Given that our interest is in the subsequent 
dynamics of the star streams the assumption is appropri¬ 
ate for the time being, but it will eventually be important 
to incorporate a more complete dynamical description of 
the originating star cluster. 

The evolution of a very low mass, collisionless, satel¬ 
lite orbiting in an external potential in which dynamical 
friction is not a factor is well handled with multipole 
expansion codes. The relatively low cost and good ac¬ 
curacy of multipole n-body codes has a rich history as 
the first codes that could handle large numbers of self- 
gravitating collisionless particles (|Aarseth||1967[ |van Al- 

J eebies|198Ut|VilJ 


bada fc van Gorkom|1977l|Fry & Peebles|1980|| Villum sen 


1982[ White 11983 McGlynn 1984). Multipole methods 

remain popular for moderately distorted spheres with ar¬ 
bitrary density profiles for whi ch there are now effi cient 
parallel processing algorithms (Meiron et al.| 2014|) and 
have bee n generalized to the pow erful self-consistent field 
method (Hernquist & Ostriker 1992). A multipole field 
provides substantiaTsmoothing of the gravitational field 
and completely suppresses the complex astrophysics of 
the two and more body interactions in a dense globu¬ 
lar cluster. If desired, star-star collisi onal p hysics can 
be added using a Monte Carlo method (Giersz 2006) al¬ 
though none is included here. 

After testing potentials calculated with up to fourth 
order multipole expansions we found that a simple shell 
code, which only calculates the spherical part of the grav¬ 
itational field, gave results that were essentially indistin¬ 
guishable, with the shell code being less subject to the 
problem of defining the center. This is not unexpected in 
the situation we are modeling since much of the mass of 
the satellite is in the central regions where the tidal field 
is weak. On the other hand, the degree of quantitative ac¬ 
curacy needs to be verified. The shell code calculates the 
radial gravitational acceleration of a particle at radius, r 
with a mildly softened gravity, —GM{< r)r/(r 2 + s 2 ) 3 / 2 
where M(< r) is simply the total mass interior to r rel¬ 
ative to a central point. The softening s is set at 20% of 
the core radius of the King model cluster we wish to sim¬ 
ulate and is needed for the case where two particles are 
near the coordinate center. Making a good choice for the 
central point of the potential is important, since a poor 
choice can lead to undesirable oscillations and even an 
unstable cluster. For the situations of interest here the 
simple approach of integrating the orbit of a fictitious 
particle placed at the position of the coordinate center 
of the cluster works extremely well and shows no signs 
of instability. 

The calculation sorts particles in order of ascending 
radii and then calculates the gravitational field at each 
successive particle, simply as the mass interior to its ra¬ 


dius. The sub-halo gravity is calculated with a direct 
sum. Our shell code has fixed time steps, At = 0.004, 
that accurately capture the tidal heating in the star clus¬ 
ter and the interactions with sub-halos down the tidal 
stream. For many of our simulations about half of the 
particles remain in the cluster so increasing the time step 
size for particles outside the cluster would not necessar¬ 
ily give rise to a large speed increase. The code was 
parallelized with a parallel quicksort and MPI routines 
to reduce the execution time. Other than the quicksort 
the algorithm is essentially trivially parallel. 

3.2. Comparison of n-body and shell code results 

As tests of the accuracy of the shell code, Figure [7] 
shows the mass as a function of time in Gadget2 and the 
shell code for a cluster started at r = 2.2 in a spherical 
isochrone. Note that difference between the full n-body 
code and the shell code is almost entirely one of slightly 
different orbital periods, not of mass lost per orbit. The 
shell code uses the high accuracy orbit of the center of 
the star cluster. The angular momentum and the radial 
action calculated at the end of the simulation when the 
potential is again spherical are measured relative to the 
star cluster as a function of their time of ejection from 
the cluster. The resulting distributions are shown for the 
n-body simulation in Figure [8] and for the shell code in 
Figure [9| The two are essentially identical, other than 
the slight difference in the orbital period, with the shell 
code having the correct period. We conclude that the 
shell code gives accurate results for the range of satellites, 
orbits and potentials of interest here. 

3.3. Streams in Triaxial Potentials without Sub-Halos 

Clusters were started at x = 4 with L/L c = 0.7 
and 0.4 in spherical and [a y ,a z \ = [—0.1,—0.05] triaxial 
isochrones with no sub-halos. A cluster in a low eccen¬ 
tricity orbit initiated with L/L c = 0.7 loses 39% of their 
mass over the interval to time 600 nearly independent 
of triaxiality. The cluster undertakes approximately 17 
radial orbits in the spherical potential. Over the same 
time a cluster on a high eccentricity orbit loses 99.7% 
of its mass in the the triaxial potential and 94% in the 
spherical potential (21 radial orbits). Although the dif¬ 
ference is easily measurable, when considered as a rate 
the differences with potential shape are small compared 
to the dependence on orbital eccentricity. The mass loss 
does depend on the location of the surface used to mea¬ 
sure contained mass. We use 0.03 units, approximately 
three times the initial tidal radius. 

Figures [10] and [IT] show the action variables relative to 
the star cluster oFthe particles in the stream as func¬ 
tion of their ejection age for two triaxial potentials. Two 
clear features are that the two tidal tails remain almost 
completely symmetric, largely a consequence of the small 
size of the cluster, which has a tidal radius of 0.007 units, 
about 0.3% of the orbital radius. Consequently, the tidal 
field is essentially symmetric across the cluster and the 
two tails are nearly the same, although moving away from 
the cluster in opposite directions. The second feature 
is that the radial action and angular momentum of the 
stars relative to the satellite vary slowly and smoothly 
over the duration which particles in the stream. The 
change of the angular momentum in the initial orbital 
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Fig. 10. — The changes of the actions of particles of a stream 
from a cluster started with L/L c = 0.7 in a smooth triaxial with 
[a y ,a z ] = [—0.1,—0.05] (top) and [—0.05,—0.1] (bottom) The ac¬ 
tions are measured at the end of the simulation in the spherical 
potential relative to the cluster center values. The differences are 
plotted as a function of the time since the particle was ejected from 
the cluster. The action variables are red for AJ r , blue for angu¬ 
lar momentum, A L and green for A J z whose values are scaled 
downward by a factor of 30 in this plot. 


plane, A J z , is largely a measure of the subsequent tilt 
of the orbits. A J z rises exponentially with time, the 
doubling time being approximately 4 radial orbits. The 
exponential growth of the orbit tipping is an important 
new feature that the triaxial potential brings into play 
for the tidal stream relative to an orbit started with the 
same velocities in a spherical or axisymmetric potential. 
As the lower panels of Figures [To| and [TT| show the same 
triaxility on different axes relative to the orbital plane 
does not show the exponential growth. Hence, orbital 
structure will depend on the details of the potential and 
the initial orbits of particles in the potential. 

4. STREAMS IN TRIAXIAL POTENTIALS WITH 
SUB-HALOS 

Two star clusters started with the same velocites in the 
a y = —0.1, a z = —0.05 halo of Figure s JTT1 and [To} but 
now with 1000 sub-halos containing 0.5% of the mass of 
the halo, is shown in projection (both at two times) onto 
the orbital plane in Figure [12] The addition of the sub¬ 
halos leads to gaps in the stream and helps to push the 
progenitor and stream onto a somewhat different orbits 
relative to the same potential without sub-halos. 

4.1. Stream Visibilities 


Fig. 11. — Same as Figure ^0] but with a cluster started with 
L/L c = 0.4. Note that the vertical scale is reduced relative to 

Fig.nrg 


A star cluster on a high eccentricity orbit is subject 
to stronger tidal fields at pericenter than a low eccen¬ 
tricity cluster with the same apocenter. The high ec¬ 
centricity orbit also takes the stream into regions where 
there are more sub-halos, therefore we expect high eccen¬ 
tricity streams to have more orbital spreading than low 
eccentricity streams. To examine this idea we compare 
the streams from clusters on L/L c = 0.4 and 0.7 orbits 
in Figure [I2| A cluster started with L/L c = 0.4 loses 
mass at more than twice the rate of a cluster started 
with L/L c = 0.7. Over the somewhat arbitrary simula¬ 
tion duration here the high eccentricity mass cluster is 
reduced to a few percent of its initial mass and it is only 
a coincidence of timing that the final moment is caught. 
Figure [12] shows the projected view of the two streams at 
two snapshot times. The upper panel shows the location 
of all particles, showing that the high eccentricity stream 
has many stars scattered well away from the stream. In 
currently available observational data, once match fil¬ 
tered to identify streams, it usually is difficult to identify 
stream locations where the local surface density is less 
than abou t a facto r of 5 in density below the peak values 
(Carlber g fc Grillmair| |2Q13). 

Figure l3|shows that the surface density distribution of 
high and low eccentricity streams have a different charac¬ 
ter, with low eccentricity streams having relatively more 
paricles at higher surface densities. The distribution is 
not strongly orbital phase dependent. The peak values 
are similar, with about 100 particles per pixel (each pixel 
is 0.01 units square) for both low and high eccentricity 
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Fig. 12.— Snapshots at times 580 and 600 of a low and high 
eccentricity stream. In the top panel the low eccentricity stream is 
the two blurry structures that come near the center, whereas the 
high eccentricity stream is the thinner structures that cross them 
and nearly connect. The contrast scale goes from 0 to 10 particles 
per 0.01 unit square pixel in the top panel and from 20 to 100 
particles per pixel in the bottom. 



Fig. 13.— The cumulative distribution, normalized to the to¬ 
tal, of surface densities in the low eccentricity stream (green and 
turquoise lines) and the high eccentricity stream at times 580 and 
600 (red and blue lines) in the [—0.1, —0.05] halo shown in Fig |12| 


stream, although we do note that the high eccentricity 
stream has somewhat higher peak densities. At a factor 
of 5 down from peak, 20 particles per pixel, about 25% of 
the low eccentricity stream is visible but only 5% of the 
high eccentricity stream. The lower panel of Figure [T2| 
shows the regions with surface densities above 20 parti¬ 


cles per pixel. More quantitatively, from Figure 13 


see that 95% of the high eccentricity stream is below the 
20% of peak surface brightness level, nearly independent 
of the orbital phase. Similar behavior is seen in the more 
general simulations of Ngan et al. (|2015b), the results 
here clarifying the combined roles of sub-halos and tri- 
axiality in allowing much of the stream to disperse to low 
densities. 

The visible segment of the Op hiuchus stream has a 
luminosity of 1.4 db: 0.6 x 10 3 L 0 (Bernard et al. 20141) 
and i s estimated to be on an e = 0.67 orbit (Sesar et al. 
2015). If the visible section of Ophiuchus has a similar 
threshold for visibility as the high eccentricity stream 
here, then the total luminosity of the entire Ophiuchus 
stream could be 20 times higher, a total of 3 x 10 4 L q 
T hese exploratory ideas mainly point to a scenario to 
understand the very short and low luminosity Ophiuchus 
stream which needs to be developed into a model to test 
its viability. 

5. STREAM ACTION VARIABLES 

Action variables are constants of the motion in a sta¬ 
tionary potential. If the actions have the same distribu¬ 
tion of values along the length of a stream, then they 
can be used along with the stream shape and any avail¬ 
able velocities to put strong constraints on the mass and 
shape of the galactic dark matter potential. The actions 
of stars ejected into a stream have pulses which vary sys¬ 
tematically around an orbit, but these details average 
together in position space a few orbital cycles further 
down the stream (Carlberg]]2015). As seen in Figures 10 
and [lT] a time variable triaxial potential induces smooth 
systematic changes in the potential which are symmetric 
relative to the satellite. Although such a smooth varia¬ 
tion will complicate modeling the constraining power of 
streams remains. 

An important issue is how dark matter sub-halos affect 
stream actions. The effects of a modest sub-halo fraction 


is shown in Figures 14 and 15 where 0.5% of the halo mass 
of the same model potentials as in Figures [To] and [To| is 


converted to 1000 equal mass sub-halos. The addition of 
sub-halos largely leaves the local character of the action 
distribution in place, a consequence of the probability of 
a sub-halo hitting the stream just as it emerges being 
fairly low. However, with increasing distance down the 
stream the cumulative action of the sub-halos is to create 
sharp offsets in the actions that can be up to a factor of 
two changes relative to the satellite values. However, 
the absolute size of the changes are typically only 1% of 
the underlying value so within a variation of about 1% 
the action variables along the stream can be considered 
approximately constant. 

6. LOW ECCENTRICITY STREAM KINEMATICS 

Low eccentricity streams have less differential shear 
and less morphological chan ge with orbital phase than 
high eccentricity streams (Carlberg 2015) making them 
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Fig. 14.— The actions of stream particles rel ative to the satellite 
center for the L/L c = 0.7 satellite shown in Fig. |12| as a function of 
ejection time. The triaxial potentials [a y ,a z \ = [— 0.05, —0.1] (top) 
and [—0.1,—0.05] (bottom) contain 1000 sub-halos with 0.5% of 
the mass. The A J z values are scaled down a factor of 30 to fit on 
the plot. 


better suited than high eccentricity streams as indicators 
of the sub-halo population of the galactic halo. Tidal 
streams started with L/L c = 0.7 in an a y = —0.1, 
a z = —0.05 halo with an increasing mass fraction in a 
varying number of sub-halos, ranging from 100 to 3000, 
along with zero for reference, are shown in Figures [16} 
The plots show the radial position of the particles rela¬ 
tive to radius of the orbit of the center of the star cluster 
at that azimuthal angle. The reference orbit is calcu¬ 
lated in the spherical isochrone potential at the end of 
the simulation. The figures show that the sub-halos are 
sufficiently numerous and massive to have significant en¬ 
counters with the streams. At a 1% sub-halo mass frac¬ 
tion the stream is becoming so badly contorted for any of 
the sub-halo numbers that 1% becomes an approximate 
upper limit to the sub-halo fraction that can be tolerated 
for thin streams, although this is consistent with the pre¬ 
dicted local fraction. With the mass entirely in 100 or 
300 sub-halos the stream encounters are rarer, but have 
a larger effect leading to large gaps. As the sub-halos 
become more numerous and less massive the stream is 
able to large retain its basic shape, with the sub-halos 
producing gaps in the stream, as is most clearly seen for 
the 3000 sub-halo case. 

The 1000 sub-halo case is repeated for a star cluster of 
mass m = 10 -7 in Figure [l7| which is ten times the clus¬ 
ter mass of Figure [16) The stream width and length scale 
with the tidal radius so are a factor of 2.15 times larger. 


Fig. 15. — Same as Figure but for satellites initiated on 
L/L c = 0.4 orbits. Note that the vertical scale is a factor of two 
larger. The cluster essentially dissolves by the end of the simulation 
so the tidal radius is reduced for more recently lost star particles. 


Two snapshots are shown to illustrate the variation with 
orbital phase. The lower panel of Figure [I7| can be com¬ 
pared to the second panel from the bottom of Figure [16} 
but note that the axes are not the same. Although the 
sub-halo effects remain constant, the figure also shows 
that sub-halo induced gaps remain visible in this slightly 
more than twice as wide stream from a ten times more 
massive progenitor cluster. Figure pT| also shows that the 
tidal “spurs” on the stream are much more evident par¬ 
ticularly near apocenter as a consequence of the orbital 
phase dependence of the mass ejection from the satellite 
being initiated at pericenter with the pa rticles d rifting 


away at apocenter as is the case for Pal 5 (Dehnen et al. 


2004). 

The radial and vertical velocities of the streams em- 
anating from the m = 10 -8 cluster of Figure 


16 


are 

shown in Figures [18] and 19 The effects of increasing 
mass fraction in sub-halosT/s/u over a range of N s h are 
clearly visible in the velocities of [18] and [T9 In all cases 
the amplitude of the perturbations grows significantly 
with increasing f s h which gives a better defined “gap” 
along the stream. The radial velocities mix together the 
velocities of the tidal features and have considerable or¬ 
bital phase dependence, whereas the velocities perpen¬ 
dicular to the orbital plane have little orbital phase de¬ 
pendence. The mean vertical velocity and its disper¬ 
sion about the mean are shown in Figure |20| Along the 
stream the dispersion ris es as the orbits spread out in 


Gaskins & Valluri 2008 


the aspherical potential (Helmi fe White 1999 Siegal- 


Ngan et al. ~ 2015b| ) . Further 
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Fig. 16.— The stream locations relative to radius of the orbit of 
the satellite, A r relative to the angular location of the satellite, cf). 
The particles from satellites on L/L c = 0.7 orbits with 100, 300, 
1000 and 3000 sub-halos. The mass in subhalos in each panel is 0 
at the bottom then 0.2%, 0.5%, and 1% (top). 

down the stream the stream has been out of the clus¬ 
ter for a longer time and has been subject to more sub¬ 
halo strikes. The mean velocity is knocked away from 
zero and there are local increases in the velocity disper¬ 
sion. on a length scale comparable to the scale size of 
the dark matter halo. The Hernquist radius of the sub- 
halos is 0.017y(/ sh /0.005)(1000/V^). For an orbital 
apocenter of 30 kpc, the sub-halo radius scales to 0.13 
kpc. Kinematic detection of sub-halo crossings requires 
enough velocity measurements per kiloparsec to be able 
to detect the expected velocity changes. The numbers 
can be worked out from simulations matched to a stream 
of interest, but it will be of order one hundred measure¬ 
ments of stream velocity per kiloparsec to be able to get 
started. For 0.5% of the halo mass in sub-halos, excur¬ 
sions around the mean velocities reach 0.02 units which 
is 5% of the circular velocity of 0.4 units. Or, if scaled to 
a circular velocity of 220 kms -1 means that values of 11 
kms -1 are expected in the largest excursions from the 
mean. For 0.2% mass in sub-halos the velocity excur¬ 
sions are proportionally smaller. The velocity dispersion 
in the thin parts of the streams, which are most read¬ 



-0.5 0 o's 1 

0 [radians] 


- 2-10 1 2 
0 [radians] 

Fig. 17.— Th e A r — cf) plot for a satellite ten times as massive 
as in Figure |l6| m = 10 -7 with no sub-halos at the bottom, then 
1000 sub-halos with 0.2%, 0.5%, and 1% (top) of the halo mass, 
for time 580 (top) and 600 (bottom) to illustrate the dependence 
on orbital phase with a time near apocenter at the top and near 
pericenter at the bottom. The bottom pan el is comparable to the 
second from the bottom panel of F igure [l6] although the scales are 
different. 

ily identified, is about 0.002 units, which scales to about 
1 kms -1 , so these excursions should be readily detectable 
given enough velocities along the stream. 

7. DISCUSSION AND CONCLUSIONS 

We have developed a framework which is able to repro¬ 
duce a large range of properties of tidal streams as are 
seen in galactic halos extracted from cosmological simu¬ 
lations. Basing the potential on the isochrone puts these 
simulations into a well studied regime and allows easy 
dynamical insight into the orbital mechanics at play. A 
central element of tidal stream analysis is that the distri¬ 
bution of star velocities in the stream can be accurately 
calculated, here with a spherical shell n-body code. Once 
free of the satellite potential, stars orbit freely, respond¬ 
ing to in the overall potential around the orbit and the 
presence of sub-halos orbiting within the host potential. 
A question of interest is how much do the orbital actions 
of stream stars change relative to their initial values as 
a result of the asphericity of the potential and any sub¬ 
halos. If the changes are large it would diminish the 
value of streams as probes of the potential shape. We 
find that for realistic local sub-halo mass fractions of 1% 
or less that the orbital actions changes are about 1-2% 
over approximately 15-20 radial orbits. The stars most 
distant from the progenitor cluster are the most affected, 
but the changes are small enough that streams can be 
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Fig. 18. — Same points as Figure p~6] except with the radial ve¬ 
locities relative to the orbital radial velocity at that angle. 


Fig. 19. — Same points as Figure [l6] except with the vertical 
velocities relative to the orbital vertical velocity at that angle. 


approximated, to this precision, as having invariant or¬ 
bital actions. The outcome is a desirable situation in 
which sub-halos should have easily measurable local ef¬ 
fects in velocity space, but the mean stream motions are 
not greatly disturbed. 

Tidal streams in spherical halo models are much thin¬ 
ner and simpler than those seen in more realistic cosmo 


logica l halos, whether there are sub-halos or not (Ngan 


et al. 2015b). The underlying cause is that aspherical 


potentials have orbits that are not present in a spheri¬ 
cal potential, opening up more ways for the stream to 
spread. Near the progenitor the stream is sufficiently 
young that the number of sub-halo encounters is low and 
the stream dynamics mainly reflect the change in tidal 
field around the orbit. Further down the stream, where 
the stream has been orbiting freely longer, the cumula¬ 
tive effect of sub-halos and potential asphericity builds 
up. In the case of a high ellipticity, e > 0.6, star cluster 
progenitor, the stream orbit spreading can be so large 
that the majority of the stream would be hard to iden¬ 
tify as a high density region on they sky. However, the 
dense segments of high eccentricity streams certainly are 
visible, even though they may represent only a few per¬ 
cent of the entire stream, leaving the rest as a ghost. It is 



0 [radians! 


Fig. 20. — The mean (blue) and dis persion (red) in vertical ve¬ 
locities along the stream of Figure [T9] for 1000 sub-halos, with the 
sub-halo mass fraction rising from 0 at the bottom to 0.2, 0.5 and 
1% at the top. 
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therefore possible that a short stream segment like Ophi- 
uchus could be the high density region of a stream with 
a total luminosity a factor of ten or more larger. 

Velocity measurements are the key to making progress 
with stream dynamics as indicators of both the shape of 
the halo and the substructure in the halo. It is now 
straightforward to reliably predict the velocity distri¬ 
bution of the particles newly injected into the stream, 
provided that the mass loss is predominantly driven by 
galactic tides. The presence of sub-halos leads to patches 
along a stream comparable in size to the characteristic 
radius of a sub-halo, typically 0.1-0.3 kpc, with an en¬ 
hanced velocity dispersion and velocity offsets from the 
mean value for the stream, as shown in Figure [20l It 
will likely be easier to detect the mean velocity offsets 
along the stream, rather than local velocity dispersion 
increases. The large scale shape of the galactic potential 
leads to smooth variations along the stream. Although 
there already is considerable incentive to acquire spec¬ 


troscopy of stars in stellar streams for abundances and 
orbital analysis, these simulations show that once large 
numbers of precision velocities are available for stream 
stars they will provide definitive tests for the presence of 
sub-halos. 

This paper has concentrated on how sub-halos and tri- 
axial potentials effect action variables, orbits and veloc¬ 
ities and left aside an examination of gaps in streams 


|Yoon, Johnston & Hogg||2011| [Carlberg 20l2l |Carlberg 

' Crillmairl2013l |ErkaI^~Beloinrov|2015| r as indicators 

of the sub-halo population. Caps will be most useful 
as indicators in low eccentricity streams, being blurred 
outjairly quickly in high eccentricity streams (Carlberg 
2015). An analysis of gaps and velocity offsets in these 


streams as sub-halo indicators will be presented else¬ 
where. 
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